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Intrinsic finite element modeling of a linear 
membrane shell problem 



Peter Hansbo * Mats G. Larson ' 



Abstract 



A Galerkin finite element method for the membrane elasticity problem on a 
meshed surface is constructed by using two-dimensional elements extended into three 
dimensions. The membrane finite element model is established using the intrinsic ap- 
proach suggested by Delfour and Zolesio [8]. 



1 Introduction 



Models of thin-shell structures are often established using differential geometry to define 
the governing differential equations in two dimensions, cf. Ciarlet [I] for an overview. A 
simpler approach is the classical engineering trick of viewing the shell as an assembly of 
flat elements, in which simple transformations of the two-dimensional stiffness matrices are 
performed, cf., e.g., Zienkiewciz [15]. In contrast to these approaches, Delfour and Zolesio 
[HI |9], [10] established elasticity models on surfaces using the signed distance function, which 

(SI can be used to describe the geometric properties of a surface. In particular, the intrinsic 

tangential derivatives were used for modeling purposes as the main differential geometric 
tool and the partial differential equations were established in three dimensions. A similar 

S^ concept had been used earlier in a finite element setting for the numerical discretization of 

the Laplace-Beltrami operator on surfaces by Dziuk |12| . resulting in a remarkably clean 
and simple implementation. For diffusion-like problems, the intrinsic approach has become 
the focal point of resent research on numerical solutions of problems posed on surfaces, cf., 

e. g ., m la rami Halm 

The purpose of this paper is to begin to explore the possibilities of the intrinsic approach 
in finite element modeling of thin-shell structures, focusing on the simplest model, that 
of the membrane shell without bending stiffness. We derive a membrane model using the 
intrinsic framework and generalize the finite element approach of [12] . Finally, we give 
some elementary numerical examples. 
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2 The membrane shell model problem 

2.1 Basic notation 

We begin by recalling the fundamentals of the approach of Delfour and Zolesio [HI El ITU] . 
Let E be a smooth two-dimensional surface imbedded in IR 3 , with outward pointing normal 
n. If we denote the signed distance function relative to S by d(x), for x e IR 3 , fulfilling 
Vd = n, we can define the domain occupied by the membrane by 

tt t = {xeR 3 : \d(x)\ <t/2}, 

where t is the thickness of the membrane. The closest point projection p : Q t — > S is given 

by 

p(x) = x — d(x)n(x), 

the Jacobian matrix of which is 

Vp = I — dV <8> n — n (g) n 

where I is the identity and ® denotes exterior product. The corresponding linear projector 
Ps = P~e(x), onto the tangent plane of £ at x G S, is given by 

P s := I - n <g> n, 

and we can then define the surface gradient Vs as 

V S :=P S V. (2.1) 

The surface gradient thus has three components, which we shall denote by 
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d 

dxf 

d 

dxf 



For a vector valued function v(x), we define the tangential Jacobian matrix as 

dvi dv\ dv\ 
dxf dxf <9x s 



v :- 



dv-2 dvo 
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dv 2 



dxf dxf dxf 
dv3 dv 3 dvs 



dxf dxf dxf 
and the surface divergence Vs • v := trVs <8> v. 



2.2 The surface strain and stress tensors 



We next define a surface strain tensor 

1 



e s («) := - (V s <8> u + (V s <g> w) T ) 



which is extensively used in [HI El HD] , where it is employed to derive models of shells based 
on purely mathematical arguments. 

From a mechanical point of view, the problem of using £s(u) as a fundamental measure 
of strain on a surface lies in it not being an in-plane tensor, in that e-^{u) n ^ 0. The shear 
strains associated with the out-of-plane direction are typically neglected in mechanical 



models, but are present in £^(u) (cf. Remark 2.1). To obtain an in-plane strain tensor we 
need to use the projection twice to define 



e£(u) 



P^e(u)Ps 



which lacks all out-of-plane strain components. For a shell, where plane stress is assumed, 
this strain tensor can still be used, since out-of-plane strains do not contribute to the strain 
energy. 



Remark 2.1 It is instructive to work out the details at a surface point whose surrounding 
is tangential to the X\X 2 -plane. In this case n = (0, 0, 1), 



and 
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The terms in e^ not present in e^ are shear strains that are typically neglected for thin 
structures, and it is clear that in our case e£ is the relevant strain tensor. 



However, the tensor e^ is rather cumbersome to use directly in a numerical implemen- 
tation; it would be much easier to work with e^ which can be establish using tangential 
derivatives. For this reason, we use the fact that there also holds (as is easily confirmed) 

e£(tt) = P S £eW-Ps = ^ ( p sV E ® MPs + (P S V E ® «P S ) T ) , 
and since n • £s(w) • n = we have the following relation: 

e£(u) = es(w) - ((es(w) • n) n + n <g> (es(w) • n)) , 
so that, using dyadic double-dot product, 



where a is a tensor and u, v are vectors, we arrive at 

e£(t*) : e£(v) = e E («) : e E (v) - 2(e E (u) • n) • (e E (v) • n), (2.2) 

which will be used in the finite element implementation below. We also note that there 
holds 

tre||(iO = VE-v, (2.3) 

where tre = J2k £ kk- 

We shall assume an isotropic stress-strain relation, 

<t = 2/xe + Atre J, 

where a is the stress tensor and J is the identity tensor. The Lame parameters A and \x 
are related to Young's modulus E and Poisson's ratio v via 

E x Eu 



2(l + i/)' (l + i/)(l-2i/)" 

For the in-plane stress tensor we thus assume 

cr£ := 2^2 + Atre^ P s , 

in the plane strain case and, in the plane stress case, which is appropriate for a thin 
membrane, 

<r£ :=2//e£ + A tre£P E , (2.4) 

where 

2A^ Eu 

Ao : 



A + 2/i 1 - u 



2 ' 



2.3 The membrane shell equations 

Consider a potential energy functional given by 



n(t*t) := - / <r{u t ) : e(u t )dtt t - f t U t 
* Jilt Jtt t 



where f t is of the form ft = f°P- Under the assumption of small thickness, we have 

-t/2 



/ /(a>)dfi t « !' I fdHdz 

J fit J -til JT, 



and thus 



U(u t ) « nf («) := f - J <r£(u) : e£(u)d£ -tff-udE 

■= jKW. e s(«))s - *(/, u) E . 

Minimizing the potential energy leads to the variational problem of finding u G V, 
where V is an appropriate Hilbert space which we specify below, such that 

a E (u,tO = te(tO V«G7 (2.5) 



where, by (2.2) and (2.3), 



a s {u,v) = (2/xe£(u),e£(u)) s + (A tre|j(u),trejj(v))E 

= (2//e s (w), e s (v)) s - (4^e E (ti) • n, e E (u) ■ n) E + (A V s • u, V E • v)e> 

and Zs(v) = (f,v)-£. This variational problem formally coincides with the one analyzed 
in the classical differential geometric setting by Ciarlet and co-workers [HI E] , as shown in 

ma. 

Splitting the displacement into a normal part u n := u ■ n and a tangential part u t := 
it — u n n we have the identity 

e£(u) = e£(u t ) + M « v ® n = e T,( u t) + m„k, 

where n = V <8> Vd is the Hessian of the distance function d, cf. [10], The bilinear form 
can therefore also be written in the form 

a s (u, v) = (2//(ef (t*t) + w„k), e£(v t ) + v n n)^ 

+ (A (tr e E (u t ) + w n tr k), tr e£(v t ) + w « tr K )s ( 2 -6) 



This means that we do not have full ellipticity in our problem. Based on this observation 
we conclude that the natural function space for the variational formulation is 



V = {v : v n e L 2 (S) and v t 6 [H\Z)} 2 }, 



cf . [5] . The loss of ellipticity have consequences for the numerics and we comment on this 
in the numerical examples below. 
Since 

(o-£(u),e£(u)) E = (er£(u),e s (u)) s 

we find, using Green's formula, the pointwise equilibrium equation 

-V E -<r£(tt) = / inS, (2.7) 



which together with the constitutive law (2.4) defines the intrinsic differential equations of 
linear elasticity on surfaces. 

3 The finite element method 

3.1 Parametrization 

Let l~h '■= {T} be a conforming, shape regular triangulation of E, resulting in a discrete 
surface £/>. We shall here consider an isoparametric parametrization of the surface (the 
same idea can however be used for arbitrary parametrizations). In the numerical examples 
below we use a piecewise linear approximation, meaning that the elements T will be planar. 
For the parametrization we wish to define a map F : (£, if) — > (x, y, z) from a reference 
triangle T defined in a local coordinate system (£, rf) to T, for all T. To this end, we 
write x = x(£,rj), where x = (x,y,z) are the physical coordinates on S^. For any given 
parametrization, we can extend it outside the surface by defining 

x(£,V,0 =B(Z,v) + tn{Z,v) 

where n is the normal and —t/2 < £ < t/2. In some models, where the surface is an 
idealized thin structure, it is natural to think of t as a thickness. 

For the representation of the geometry, we first introduce the following approximation 

of the normal: 

n h 



n~n h = p£r> n o = Yl n ^(^ V)* 



\n 



where <Pi(£,7)) are the finite element shape functions on the reference element (assumed 
linear in this paper), and «,j denotes the normals in the nodes of the mesh. We then 
consider parametrizations of the type 

x(£, v, C) » x h (t, V,0 = Y, fawfo V) + C nupifo v)) (3.1) 

i 

where X{ are the physical location of the nodes on the surface. For the approximation of 
the solution, we use a constant extension, 

uwu'^^il) (3-2) 



where Ui are the nodal displacements, so that the finite element method is, in a sense, 
superparametric. Note that only the in-plane variation of the approximate solution will 
matter since we are looking at in-plane stresses and strains. We employ the usual finite 
element approximation of the physical derivatives of the chosen basis {fi} on the surface, 
at (£,//), as 
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d( Pj 

§ y 

. dz . 

This gives, at ( 
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With the approximate normals we explicitly obtain 
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Remark 3.1 The approach by Dziuk [12] (and also the classical engineering approach, 
fWj) is, in our setting, a constant- by- element extension of the geometry using facet triangles 



{T} so that, with n? the normal to the facet, x (£,?7,C)| T 

dip., 



d(fj 
dx 
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This low order approximation has the advantage of yielding a constant Jacobian from a 
linear approximation. For some applications this is, however, offset by the problem of 
having a discontinuous normal between elements. 



3.2 Finite element formulation 

We can now introduce finite element spaces constructed from the basis previously discussed 
by defining 

W h := {v:v\ T oFe P k (f), VT G T h ; v G C°(£ ft )}, (3.3) 

(in the numerics, we use k — 1), and the finite element method reads: Find Uh G V h := 
[W h f such that 

av h {u h ,v) = l* h {v), VveV\ (3.4) 

where 

ax h (u, v) = (2jue Eh (w), £s h (v))s h - (4jue E „(u) • n 7 *, e E/ » • n h ) E „ 

+ (A V s ,-w,V Sh -v) Sh 

and/ Sfc (v) = (f,v) Sh . 

3.3 Extension to surfaces with a boundary 

If the surface £ has a boundary <9£ we assume that <9E = Uj<9£j where <9£j are closed 
components. On each of the components <9Ej we let qr : <9Xj — >■ R 3 , j = 1, 2, 3, be smooth 
orthonormal vector fields. We strongly impose homogeneous Dirichlet boundary conditions 
of the type 

qj ■ u = on <9Ej, 1 < j < di, (3.5) 

where di = 1, 2, or 3, and weakly the remaining Neumann condition 

(n d x ■ <r£(iO) • <?j = 0, di< j < 3, (3.6) 

where n^s is the unit vector that is normal to <9£ and tangent to E. Note that not every 
combination of boundary conditions and right hand side leads to a well posed problem. 

4 Numerical examples 

In the numerical examples below, the geometry is represented by flat facets, and the 
normals are taken as the exact normal in the nodes, interpolated linearly inside each 
element. Our experience is that similar results are obtained if we use L 2 — projections of 
the flat facet normals in the nodes and then interpolate these linearly. 

4.1 Pulling a cylinder 

We consider a cylindrical shell of radius r and thickness t, with open ends at x = and at 
x = L, and with fixed longitudinal displacements at x — 0, and radial at x — L, carrying 
a horizontal surface load per unit area 

f( \ F x 

Zirr L z 



where F has the unit of force. The resulting longitudinal stress is 

F (1 - {x/Lf) 



a 



Anrt 



We take as an example a cylinder of radius r = 1 and length L — 4, with material 
data E = 100 and v = 1/2, with thickness t = 10~ 2 , and with F = 1. In Fig. [I] we 
show the solution (exaggerated 10 times) on a particular mesh (shown in Fig. [2]). Note 
that the lateral contraction creates radial displacements depending on the size of stress. 
Finally, in Fig. p we show the L 2 error in stresses, \\cr — crh\\L 2 (n), where <r := cr^iu) and 
(Th := cr^(uh), which shows the expected first order convergence for our P 1 approximation. 
The black triangle shows the 1:1 slope. 



4.2 A torus with internal pressure 

We consider a torus with internal gauge pressure p for which the stresses are statically 
determinate. Using the angle and radii defined in Fig. |4j the principal stresses are given 
by 

pr pr ( r sin ( 

°"l = 7T7; °"2 = — 1 - 



2t t V 2(R + rsin6 

where o\ is the longitudinal stress, a 2 the hoop stress, and t is the thickness of the surface 
of the torus. The constitutive parameters and thickness where chosen as in the cylinder 
example, and we set R = 1, r = 1/2, and p — 1. 

Again we compute the stress error ||er — 0"h||z a (fi). We show the observed convergence 
in Fig. [7] at a rate of about 3/4 (the slope of the black triangle), which is suboptimal, but 
does occur in problems where elliptic regularity is an issue, cf. [3j, Lemma 10. We thus 
attribute this loss of convergence to the load now being in the normal direction of the shell, 
for which we do not have ellipticity. 
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Figure 1: Displacements (exaggerated by one order of magnitude) on a particular mesh. 
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Figure 2: The cylinder before deformation. 
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Figure 3: Stress convergence for the cylinder 
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Figure 4: A cut through of the torus 




Figure 5: A typical mesh on the torus. 
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Figure 6: Deformations in the torus case, exaggerated by two orders of magnitude. 
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Figure 7: Convergence of the stresses in the torus case. 
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